//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
/** @author  John A. Miller
 *  @author  Michael E. Cotterell
 *  @see     LICENSE (MIT style license file).
 */

package simopt

import collection.mutable.ArrayBuffer

import scalation.linalgebra.{ MatrixD, VectorD }
import scalation.linalgebra_gen.Vectors.VectorI
import scalation.metamodel.QuadraticFit
import scalation.optima.IntegerTabuSearch

class ERRegressionModel (u: (VectorD) => Double) extends ERParams
{
    val qf = new QuadraticFit (u, 5) // construct a quadratic fit object 
    var current: VectorI = null      // the current input vector
    var grid = null                  // the design grid

    // returns true if x is in the current quadrant
    def sameQuadrant (x: VectorI): Boolean = {
        if (current == null) false
        else (0 until x.dim) map { i => i match {
                case _ if (x(i) <= 5 && current(i) <= 5) => true
                case _ if (x(i) >  5 && current(i) >  5) => true
                case _ => false
            } // match
        } reduceLeft (_ && _)
    } // sameQuadrant

    // returns the design grid for the quadrant containing x
    def designGrid (x: VectorI): ArrayBuffer [VectorD] = {

        var grid  = ArrayBuffer.empty[VectorD]  // the grid
        var parts = Array.ofDim[VectorI](x.dim) // the parts of each dimension

        // each dimension has 3 points depending on where x is
        for (i <- 0 until x.dim) {
            if (x(i) <= 5) parts(i) = new VectorI (1, 3, 5)
            else parts(i) = new VectorI (7, 9, 11) 
        } // for

        // recursively accumulates vec while exhausting rest, adding vec to grid
        // once rest is exhausted 
        def populate (rest: Array[VectorI] = parts, vec: VectorI = new VectorI (0)) {
            if (rest.length == 0) grid += vec
            else for (point <- rest(0)) populate (rest.slice(1, rest.length), vec ++ point)
        } // populate

        populate () // populate the grid
        grid        // return the design grid

    } // designGrid

    def eval (x: VectorI): Double = {

        for (i <- 0 until x.dim if x(i) < 1)  return Double.NegativeInfinity
        for (i <- 0 until x.dim if x(i) > 10) return Double.NegativeInfinity

        // if we're in a new quadrant, reform and refit
        if (!sameQuadrant(x)) { 
            val design = qf.response(designGrid(x)) // design matrix and response vector
            qf.fit (design)                         // fit surface
        } // if

        current = x                // update current coordinate
        qf.qFormsEval(x.toDouble)  // estimate a value for the function at x

    } // u_r

} // ERRegressionModel

object ERRegressionTest1 extends App with EROpt
{
    def u_r (x: VectorD): Double = _u_p (x.toVectorI)       // use process interaction
    val m           = new ERRegressionModel (u_r)           // regression model
    val x0          = VectorI (1, 1, 1, 1, 1)               // starting point
    val its         = new IntegerTabuSearch (m.eval)        // setup optimizer
    val (x, result) = its.maximize (x0)                     // get results
    println ("x = %s; m.eval value = %f".format(x, result)) // print results
} // ERRegressionTest1

object ERRegressionTest2 extends App with ERQueueingOpt
    def u_r (x: VectorD): Double = _u_q (x.toVectorI)       // use queueing model
    val m           = new ERRegressionModel (u_r)           // regression model
    val x0          = VectorI (1, 1, 1, 1, 1)               // starting point
    val its         = new IntegerTabuSearch (m.eval)        // setup optimizer
    val (x, result) = its.maximize (x0)                     // get results
    println ("x = %s; m.eval value = %f".format(x, result)) // print results
} // ERegressionTest2
